This section is intended to introduce various aspects of the art, which may be associated with exemplary embodiments of the present techniques. This discussion is believed to assist in providing a framework to facilitate a better understanding of particular aspects of the present techniques. Accordingly, it should be understood that this section should be read in this light, and not necessarily as admissions of prior art.
Modern society is greatly dependant on the use of hydrocarbons for fuels and chemical feedstocks. Hydrocarbons are generally found in subsurface rock formations that can be termed “reservoirs.” Removing hydrocarbons from the reservoirs depends on numerous physical properties of the rock formations, such as the permeability of the rock containing the hydrocarbons, the ability of the hydrocarbons to flow through the rock formations, and the proportion of hydrocarbons present, among others.
Often, mathematical models termed “simulation models” are used to simulate hydrocarbon reservoirs and optimize the production of the hydrocarbons. A simulation model is a type of computational fluid dynamics simulation where a set of partial differential equations (PDE's) which govern multi-phase, multi-component fluid flow through porous media and the connected facility network is approximated and solved. This is an iterative, time-stepping process where a particular hydrocarbon production strategy is optimized.
Simulation models map the underlying partial differential equations on a structured (or unstructured) grid, which may be termed discretization. The grid, which may be termed a computational mesh, represents the reservoir rock, wells, and surface facility network. State variables, such as pressure and saturation, are defined at each grid block. The goal of a simulation model is generally to understand the flow patterns of the underlying geology in order to optimize the production of hydrocarbons from a set of wells and surface facilities. During the past five decades, the size and complexity of simulation models have grown proportionally with the increased availability of computing capacity. Complex simulation models often require the use of parallel computing systems and algorithms to provide adequate simulation turnaround time.
Reservoir simulation is of great interest because it infers the behavior of a real hydrocarbon-bearing reservoir from the performance of a model of that reservoir. The typical objective of reservoir simulation is to understand the complex chemical, physical and fluid flow processes occurring in the reservoir sufficiently well to predict future behavior of the reservoir to maximize hydrocarbon recovery. Reservoir simulation may often refer to the hydrodynamics of flow within a reservoir, but, in a larger sense, reservoir simulation can also refer to a total hydrocarbon system. This may include not only the reservoir, but also injection and/or production wells, surface flow lines, and surface processing facilities. Reservoir simulation calculations in such hydrocarbon systems can be based on simulations of fluid flow throughout the entire system. These calculations are performed with varying degrees of rigor, depending on the requirements of the particular simulation study and the capabilities of the simulation software being used.
The principle of numerical simulation is to numerically solve equations describing a physical phenomenon using a computer. Such equations are generally algebraic equations, ordinary differential equations (ODE), and/or partial differential equations (PDE). As a means for solving differential equations, ODE and/or PDE, numerically, there are known methods such that the finite difference method, the finite volume method, the finite element method, and the like. Regardless of which method is used, the physical system to be modeled is divided into cells, a set of which is called a grid or mesh, and the state variables that vary in time and in space throughout the model are represented by sets of values for each cell. The reservoir rock properties such as porosity and permeability are typically assumed constant inside a cell. Other variables such as fluid pressure and phase saturation are defined at specified points, sometime called nodes, within the cell (or sometimes on the boundaries between the cells). A link between two nodes is called a “connection.” Fluid flow between two cells is typically modeled as flow along the connection between them.
Most reservoir simulators use so-called structured grids in which the cells are assumed to be three-dimensional rectangular shapes distorted to conform as well as possible to geological features and flow patterns. Certain geological features and modeling situations cannot be represented well by structured grids. For example, slanted wells, faults, or fractures, which can naturally occur in the reservoir models, are almost impossible to model using structured grids.
These shortcomings can be overcome in part by using local refinement, in which selected cells are subdivided into smaller cells, and non-neighbor connections, which allow flow between cells that are physically adjacent to each other but are not adjacent in the data structure. A more powerful solution to this problem is to exploit the flexibility provided by layered unstructured grids, when a computational domain is split into geological layers and, in each layer, a grid is formed by a number of laterally contiguous grid cells. The cells forming any layer have corresponding neighbors in the layers located above and below the current one. Such a grid usually is called “2.5D grid.”
Even such grids are not sufficient to describe all the variability of complex geologic structures present in real hydrocarbon reservoirs such as branching wells, complex Y-faults, or pinch-outs, i.e. cases when the geological layer disappear inside the reservoir. The only solution to this problem is to exploit the flexibility provided by an unstructured grid.
In a simulation, a set of equations is developed to express the fundamental principles of conservation of mass, energy, and/or momentum within each cell and of movement of mass, energy, and/or momentum between cells. The replacement of the state variables that vary in space throughout a model by a finite number of variable values for each cell is called “discretization.” For proper resolution of the modeled physical processes, hundreds of thousands and even millions of cells may be required. That may lead to the millions of equations to be solved.
In order to analyze a phenomenon changing in time, it is necessary to calculate physical quantities at discrete intervals of time called time steps, irrespective of the continuously changing conditions as a function of time. Time-dependent modeling of the transport processes therefore proceeds in a sequence of time steps. During a time step, transport of various kinds occurs between cells. Through this transport, a cell can exchange mass, momentum, or energy with other nearby cells.
The equations governing the behavior of each cell during a time step couple the mass, momentum, and energy conservation principles to the transport equations. At every time step, the simulator must solve one, or more, large matrix equations, Ax=b, where A is a matrix, b is a known right hand side vector, and x is a vector of unknowns, with the number of unknowns depending on the type of time step computation method being used. Because matrix equations are quite large, having at least one equation per cell, they are solved iteratively except in case of small models.
Various time step computations have been proposed for reservoir simulation. Two most commonly used calculation methods are called “IMPES” and “Fully Implicit.” In the IMPES method, which is derived from the term “implicit-pressure, explicit-saturation,” flows between neighboring cells are computed based on pressures at their values at the end of the time step. The pressures at the end of the IMPES time step are interdependent and must be determined simultaneously. This method is called “implicit” because each pressure depends on other quantities that are known only implicitly. The basic procedure is to form a matrix equation, Axp=bp, that is implicit in pressures only, solve the matrix equation for the pressures xp, and then use the pressures in computing saturations explicitly cell by cell. In this fashion, after the pressures have been advanced in time, the saturations are updated explicitly. After the saturations are calculated, updated physical properties such as relative permeability's and capillary pressures can be calculated. Those are explicitly used at the next time step. Similar treatment can be used for other possible solution variables such as concentrations, component masses, temperature, or internal energy.
The fully implicit method treats both pressure and saturations implicitly. Flow rates are computed using phase pressures and saturations at the end of each time step. The calculation of flow rates, pressures, and saturations involves the solution of nonlinear equations using a suitable iterative technique. At each iteration, the method constructs and solves a matrix equation, the unknowns of which (pressure and saturations) change over the iteration:
            [                                                  A              pp                                                          A              ps                                                                          A              sp                                                          A              ss                                          ]        ⁡          [                                                  x              p                                                                          x              s                                          ]        =            [                                                  b              p                                                                          b              s                                          ]        .  
The matrix of this matrix equation is often called “Jacobian.” As the pressures and saturations are solved, the updating of these terms continues using new values of pressure and saturation. The iteration process terminates when predetermined convergence criteria are satisfied. Further information about reservoir simulation and computation techniques may be found in D. W. Peaceman, Fundamentals of Numerical Reservoir Simulation, Elsevier, N.Y., 1977; K. Aziz and A. Settari, Petroleum Reservoir Simulation, Applied Science Publishers Ltd, Barking, Essex, England, 1979; C. C. Mattaxand R. L. Dalton, Reservoir Simulation, Monograph, Volume 13, Society of Petroleum Engineers, 1990; U.S. Pat. No. 6,052,520 to J. W. Watts; and U.S. Pat. No. 6,662,146 to J. W. Watts. As mentioned above the matrix equations obtained after discretization of the governing equations for reservoir simulation may have hundreds of thousands and even millions of unknowns. In practice, such systems can be solved only by iterative methods. See, for example, Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd ed., Society for Industrial and Applied Mathematics (SIAM), 2003. The convergence rate of iterative methods and hence, the time spent on solving such systems, depend on the property of the matrices which is called condition number. The usual way to reduce the computational time is to reduce the condition number by employing a preconditioner B that is rewriting the system using the equation:B−1Ax=B−1b. In this equation, matrix B provides reasonable approximations to matrix A. Further, matrix B can be inexpensively inverted.
There are many different preconditioners for preconditioned iterative methods to solve matrix equations arising in reservoir simulation. Some of them can be applied only to the problems approximated on structured (See, e.g., P. T. Woo, S. J. Roberts, and F. G. Gustavson, Application of Sparse Matrix Techniques in Reservoir Simulation, SPE 4544, 1973) or almost structured grids (European Patent Application EP 2 068 263) and hence cannot be applied in case of generic unstructured grids. Some other methods apply incomplete LU factorizations (U.S. Pat. No. 6,230,101 B1 to J. R. Wallis) or singular value decomposition (SVD) reduction (U.S. Patent Application Publication No. 2007/0255779), both hardly parallelizable. Algebraic multigrid (AMG) methods are known to be most efficient methods of solving reservoir matrix equations. T. Clees and L. Ganzer, An Efficient Algebraic Multigrid Solver Strategy for Adaptive Implicit Methods in Oil Reservoir Simulation, SPE 105789, 2007. Unfortunately, the coarsening process of the AMG method applied directly to the reservoir matrix A may lead to the matrices with the increased nonzero patterns on the coarser levels, especially for non-symmetric problems, which can be not very useful for efficient parallelization. Much better approach is to combine multiplicatively a version of block-ILU0 preconditioner for the full system with an AMG cycle for the pressure block. T. M. Al-Shaalan, H. Klie, A. H. Dogru, and M. F. Wheeler, Studies of Robust Two Stage Preconditioners for the Solution of Fully Implicit Multiphase Flow Problems, SPE 118722, 2009; J. M. Gratien, T. Guignon, J. F. Magras, P. Quandalle, and O. Ricois, Scalability and Load-Balancing Problems in Parallel Reservoir Simulation, SPE 106023, 2007. This combinative approach is known as the “constrained pressure residual” (CPR) technique.
Because general-purpose computing in reservoir simulation can require substantial computing resources, proposals have been made to subdivide a simulation model into smaller segments and to perform computations in parallel on multi-processor computer or cluster of (possibly multi-processor) computing nodes. The principal attraction of parallel computing is the ability to reduce the elapsed time of simulation, ideally by a factor of N for an N-processor computing system. Parallel computing falls short of the ideal because of several factors, including recursion in linear equation solution, the overhead associated with message passing required for various computations, and load imbalances due to heterogeneities in the problem physics and characterization of the hydrocarbon fluids.
There are several approaches to constructing parallel preconditioners for linear problems arising in reservoir simulation. Most of them are based on overlapping domain decomposition. O. Dubois, I. D. Mishev, and L. Zikatanov, Energy Minimizing Bases for Efficient Multiscale Modeling and Linear Solvers in Reservoir Simulation, SPE 119137, 2009; U.S. Pat. No. 7,516,056 to J. Wallis, et al.; P. Lu, J. S. Shaw, T. K. Eccles, I. D. Mishev, and B. L. Beckner, Experience with Numerical Stability, Formulation and Parallel Efficiency of Adaptive Implicit Methods, SPE 118977, 2009. The basic idea of this method is that the computational domain is split into (non-overlapping) subdomains and each subdomain is assigned to a different processing unit of a parallel computer. Then each subdomain is extended by one or several layers of neighboring cells to get overlap. On each subdomain, the reservoir matrix is restricted to the cells of that (extended) subdomain. The local problems then can be computed in parallel by different processors using any existing technique (e.g. ILU or AMG) and the global solution is obtained by special type superposition of the local solutions from each subdomain. Generally, the convergence of the iterative method with overlapping domain decomposition is independent of the mesh size when the width of the overlapping region is kept proportional to the diameter of the subdomains. However, the convergence deteriorates linearly as the number of subdomains increases, preventing the method to be scalable. To remedy this problem, a coarse component can be added to the preconditioner, which allows for a global mean computation and communication across all subdomains. Such a two-level preconditioner is closely related to a multigrid (a two-grid) method, but in which the coarse space is generally much coarser compared to the fine grid. The performance of this type of preconditioner depends strongly on the decomposition of the domain.
U.S. Patent Application Publication No. 2007/0010979 by Wallis, et al., discloses an overlapping domain decomposition with one layer of overlap but special procedure for implementation of multigrid method. The connections between any two adjacent subdomains are lamped together with no regard of their respective values. Even though the method claims the efficiency of the obtained in such a way preconditioner is good, there is an overhead connected with small connections (by value) as they are transferred to the coarser levels.
The approach disclosed in U.S. Patent Application Publication No. 2006/0245667 by Fung, et al., can be applied to structured grid, 2.5D grid, or the combination of both as long as the resulting grid has layered structure. The method is claimed to work on a variety of computer platforms, such as shared-memory computers, distributed memory computers, or personal computer clusters. The approach heavily depends on the assumption of layered structure of the underlying grid. It splits computational domain laterally into vertical columns and distributes the split parts into separate computational nodes. The method described in that patent cannot be applied to the reservoir simulation models defined on general unstructured grids or can use any other type of partitioning of the grid cells between the computational nodes.
All the approaches described above employ a Krylov-type iterative solver like gradient minimal residual (GMRES) or a biconjugate gradient stabilized (BiCGStab). See Saad, Y., Iterative Methods for Sparse Linear Systems, 2nd ed., Society for Industrial and Applied Mathematics (SIAM), 2003. One of the properties of these solvers is that on each iteration there is a need to compute two or more inner products. In distributed parallel environment that means that for each iteration a small amount of information should be distributed globally between all the processors. However, global message passing has bad scalability, e.g., its cost grows as the number of processing units involved in the computations increases. See Message Passing Interface Forum, MPI: A Message-Passing Interface Standard, Version 2.1 (June 2008). Moreover, this global transfer of information enforces natural barriers, which may also create additional misbalance in computations.
Thus, when domain decomposition or multigrid method is used on the parallel system with multiple CPUs, the information should be exchanged on each iteration between neighbor CPUs and globally, between all CPUs. In addition, multigrid methods are based on reduction in the number of unknowns in matrix equation, which usually leads to the increased number of connections between unknowns on the coarser levels, which in turn increases the cost of exchange of the information between CPUs during each iteration.
Addressing the aforementioned issues may lead to more efficient methods of solving sparse matrix equations for reservoir models on parallel computers with multiple CPUs. Namely, an iterative method, which limits or avoids global communications between CPUs may be more efficient than methods that propagate information globally. A multigrid method that is capable of controlling the number of connections between unknowns on coarser levels will be more efficient than the one in which the number of connections grows uncontrollably.